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Abstract. We consider the dynamics of the system of self propelling particles modeled via the Vicsek algorithm in continuum 
time limit. It is shown that the alignment process for the velocities can be subdivided into two regimes: "fast" kinetic and 
"slow" hydrodynamic ones. In fast kinetic regime the alignment of the particle velocity to the local neighborhood takes 
place with characteristic relaxation time. So that the bigger regions arise with the velocity alignment. These regions align their 
velocities thus giving rise to hydrodynamic regime of the dynamics. We propose the mean-field like approach in which we take 
into account the correlations between density and velocity. The comparison of the theoretical predictions with the numerical 
simulations is given. The relation between Vicsek model in the zero velocity limit and the Kuramoto model is stated. The 
mean-field approach accounting for the dynamic change of the neighborhood is proposed. The nature of the discontinuity of 
the dependence of the order parameter in case of vectorial noise revealed in Gregorie and Chaite, Phys. Rev. Lett., 92, 025702 
(2004) is discussed and the explanation of it is proposed. 
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INTRODUCTION 

The equilibrium statistical mechanics and thermodynamics of Hamiltonian systems are well developed areas of 
Statistical Physics. There are also a lot of remarkable results for open systems, which considered far from equilibrium 
[1, 2]. In general one can consider the systems with some constraints which bounds the coordinates of the particles 
since the general formalism remains unchanged (e.g. Liouville equation etc.). It is well known that due to instability 
of trajectories the Hamiltonian systems reach the thermodynamic equilibrium which can be characterized several 
macroscopic parameters of state despite huge number of microscopic degrees of freedom. In addition total momentum 
and angular momentum are conserved. Though for every specific configuration the formation of stationary local 
vortical structures [3, 4, 5] may occur due to conservation of the angular momentum. Obviously the inclusion of 
potential forces has little grounds for the systems of intellectual particles (individuals in flocks, crowds etc.), which 
differ very much in this respect from the molecular system for which either all forces have mainly potential character 
or dissipative. 

In [6] the minimal model, the so called Vicsek model (VM) of such type was introduced. The dynamic rule for the 
alignment of the particle' velocities constructed in such a way that at high density the kinetic energy of disordered 
motion is transformed into the one of ordered motion so that the total kinetic energy is conserved. Then the system 
reaches the final state with nonzeroth total momentum even in the low (one- and two-) dimensional cases. In such a 
case the appearance of the ordered state is predetermined by the dynamic rule. Note that the VM does not take into 
account the potential interparticle forces. 

In this paper we consider the kinetic regime for the VM when the particle aligns along the velocity direction of 
its neighborhood and give the estimation for the critical noise amplitude of order-disorder transition. The structure of 
paper is as follows. In Section 1 we derive the continuum time analog of the VM and show that angular velocity of 
a particle consists of two terms which describe alignment. One of the terms describes the fast kinetic relaxation to 
the local direction the other one describes the hydrodynamic regime of alignment between the domains where local 
alignment is settled down. In Section 2 the process of the relaxation of one-particle velocity to the local value of the 
neighborhood is considered. The dependence of the rate of the relaxation on the density is obtained by the numerical 
experiment and the corresponding Fokker-Plank equation is derived. In Section 3 the influence of the dynamic nature 
of the environment is discussed. The results obtained are summarized in the conclusion. 



1. VICSEK MODEL IN CONTINUUM TIME LIMIT AND THE TWO REGIMES OF 

THE DYNAMIC 



The Vicsek model of the dynamic of self-propelling particles [6] can be represented by the relation: 



Vj(n + 1) xiij(n) =0, Vi,n, 
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is the unit vector corresponding to the averaged velocity of the neighborhood and H(r,j) is the averaging kernel with 
the characteristic averaging scale R. The absolute value of the velocity of each particle is assumed to be constant, i.e. 



| Vt(n+l) | = | M{n) |= vi 



(3) 



The noise is not included. For the VM, H is proportional to a Heaviside step function. One can also consider other 
models for the averaging kernel. One can say that the dynamics of individual particle is subjected to reduction of the 
difference between the direction of its velocity and that of the average velocity of the surrounding given by Eq. (2). 

Note that at every step the direction of the velocity v(n + 1) coincides exactly with the direction of u(n). Another 
words, the vector u(n) remains unchanged during the velocities updating. This is specific for discrete formulation but 
this is not the case in the continuous time limit since both vectors are rotating during infinitesimal time interval 8t . In 
such a limit the angular velocity o v consists of two parts. The first one is the angular velocity co u of the unit vector 
u for the average velocity of the nearest neighbors. The second one is the relative angular velocity. When the time is 
continuous taking into account constraint Eq. (1) the equation of motion of such a particle can be written as: 
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Here co v , is the "angular velocity" of ;'-th particle. 

This angular velocity depends on the velocities of neighboring particles. The self-propelling force and the frictional 
force are assumed to balance each other. The hydrodynamic model which is based on the equations of motion (1) and 
it continual analog (4) was proposed in [? ? ]. 

Note that in the discrete CVA at every step the direction of the velocity v(n + 1 ) coincides exactly with the direction 
of u(n). Another words, the vector u(«) remains unchanged during the velocities updating. This is specific for discrete 
formulation but this is not the case in the continuous time limit since both vectors are rotating during infinitesimal time 
interval 8t . In such a limit the angular velocity © v consists of two parts. The first one is the angular velocity ©„ of the 
unit vector u for the average velocity of the nearest neighbors. The second one is the relative angular velocity. 



©v, = (0 Ui + COvu , 



(5) 



where 



co Ul =U; x u, 

<B VU =AVjX Uj. 



(6) 
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The quantity A > is inverse to characteristic length scale. The latter is the radius of interaction R and is the parameter 
of the model. Indeed, in the limit R — > °o each particle has the same neighborhood, provided that N 3> 1, i.e. u, does 
not depend on ;'. Therefore, in such a limit they all has the same angular velocity, which is given by the first term of 
Eq. (5). From Eq. (5) and Eqs. (6), (7) for 2D case we obtain: 



© Vi = -A (V/ • Ui) (fl)y. - fl) u .) - (©„,. + ) . 



(8) 



2. ONE-PARTICLE RELAXATION 



Here we show that the equation (5) in 2D is closely related to well known Kuramoto model (KM) for the phase 
synchronization. Indeed, let the angle 0, characterize the direction of the velocity of ;'-th particle, then <B Vi = 0, and the 
equation (5) takes the form: 

d i = \jf i +Asm(\if i -d i ) (9) 

Here i^i denotes the angle which determines the direction vector u;. It is one of the variant of the short-range version 
of the KM [? ] (see also [10] and reference therein) of the form: 

e, = (0,+K £ sin (0,-0/) , (10) 

(i,J) 

where the brackets stand for nearest-neighbor oscillators. Thus we can state that in the zero velocity limit the Vicsek 
model with continuum time belongs to the short-range Kuramoto model class [10]. This allow to conclude that for low 
velocity the ordering in the Vicsek model is governed by the same mechanism as the synchronization in the KM. Since 
the synchronization transition is of continuum type one can expect that the continuum character of the transition take 
place for low enough velocity in Vicsek model too. This conclusion is in correspondence with the results of [1 1]. 

According to its definition vector u, changes slower than the velocity of a particle v,-. From Eq. (8) it follows 
that the second term in Eq. (5) governs the kinetic of the alignment process while the first term is of hydrodynamic 
character since it determines the behavior on scales greater than R. Therefore Eq. (8) shows that in continuous time 
limit the CVA system has the stable state where the particles align along some direction, which is characterized by the 
director uo- Equations (4), (5) can serve as the continuum time analog for the CVA. Additional confirmation of that 
is the behavior of these terms with respect to reverting the time t — > —t. The first term changes its sign and therefore 
produces the reversible contribution to the equation of motion, while the second term does not change the sign thus 
representing irreversible part of the CVA, which governs irreversible one-particle kinetics of the alignment. To study 
Eq. (8) analytically we use the approximation which takes into account that the variable u, is the "collective" one, thus 
there is the time interval which we call "kinetic" regime where it changes much slower so that co Uj and its derivative 
can be omitted. In addition we assume that the value of A does not depend on time which reflects that the number of 
"interacting" neighbors remains constant. In such an approximation Eq. (8) reduces to more simple form: 

© V; = -A(vj-Uj.) © v .. (11) 

In scalar form for the angle of the alignment a, between the vectors v, and u, taking into account that » v , = « and 

V; • U, = COS a , 

after integrating Eq. (1 1) we obtain: 

a = -A sin a. (12) 

where A is some parameter which determines the alignment rate and obviously depends on the density and the average 
velocity of the neighbors. Here we put the following initial condition a = 0, which is in accordance with that in 
simulations. The solution of Eq. (12) is: 

tan^- = tan "^" exp( — At ) , (13) 

Thus the one-particle alignment process has relaxation character. 

To compare this result with the simulation data we performed the simulation with the small initial disalignment of 

the directions of the particles in dense system p = N ( j ) w 1 . The results of simulation are represented on Fig. 1 
and demonstrate the existence of such an interval, where the dependence given by Eq. (11) takes place. 

The kinetic regime of the system subjected to the stochastic increment of the direction [6, 12] can be described by 
stochastic modification of Eq. (12): 

a = —Asiaa+L(t), (14) 

where L(t) is the standard white noise term. Then Eq. (14) is equivalent to Fokker-Plank equation for the density 
distribution function f y (a,t) [13]: 

_ = _ (Asina/v)+D __ (is) 



where D is the diffusion coefficient and we use the approximation: 

A = Xu 



(16) 



for the alignment rate for the dimension reasons, where A is the local density. Such density dependence of the 
relaxation rate is supported by the simulations (see Fig. 2). 

The case D = corresponds to the deterministic case, with the solution: 



/f (o,0 = 



tan f ; 
sin a 



(17) 



which demonstrates the alignment process in accordance with Eqs. (12), (13). The stationary solution of Eq. (15) is: 



fi st) {a)=C{D)e X T>™a 



(18) 



Note that the distribution (18) was used in [14] for the lattice model as the analog of the Boltzman distribution. Thus 
the consideration presented above can serve as the ground for such representation. 
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FIGURE 1. The log plot for the value tan ^ as the function of time (number of steps) at different radius of interaction r. The 
results of simulation with N = 900 , L = 30 , v = 0. 1 . 



3. THE ACCOUNT OF DYNAMICAL ENVIRONMENT 

In the zero velocity limit the Vicsek model can be considered in terms of the lattice model as the systems of interacting 
spins with the interaction favoring the alignment and the emergence of the long range order. Yet there is the major 
difference between the CVA and the equilibrium models. This is the coupling between the density and the velocity 
fields. Due to such coupling in the static (v = 0) case the the equilibrium systems does not order for densities below 
some threshold value close to 1 (which corresponds to the percolation threshold of randomly distributed spheres) 
while in the SSP ordering is found for all velocities [11]. The instant change of the environment in the neighborhood 
of the particle can be considered as the noise factor correlated with the local value of the order parameter - the average 
velocity of the neighbors: 
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FIGURE 2. The dependence of the relaxation rate on the number of the neighbors X = pr 2 . 



where (...) stands for the nearest neighbors of ;'-th particle and Nj is the number of such neighbors. Thus the vector u, 
is just the weighted sum of the random vectors. The number of summands is also random and describes the dynamic 
coupling between the density and the velocity of the neighbors. Let one-particle distribution function is / v . Since 
the order parameter u is the collective variable as has been said above it changes more slowly than the one-particle 
function. Thus one can consider u as the parameter for distribution function / v . Using the standard procedure analogous 
to the mean-field approach it is possible to get the self-consistent equation for the order parameter. Indeed, assuming 
that all the correlation between the particles are incorporated into u we can find the relation between the characteristic 
functions for v and u: 

G u (k) = ( e iku )=f>„G v (-Y . (20) 

Here p n is the probability of n neighbors to appear within the interaction range R and 

Gv = {e^) (21) 

is the characteristic function of the velocity distribution of a particle / v . It depends on the average value of the local 
order parameter (19). In order to get the specific results for the order parameter one need to specify the form of the 
density distribution for the number of neighbors and the distribution function / v for the velocity. The latter depends 
on the type of noise introduced into the motion. The original variant of the CVA [6, 15] used so called scalar noise 
when the direction of a particle is updated with the random increment of the angle. In recent work [16] another type 
of noise, so called vectorial noise, was proposed. We will consider these type of noise below. 

To simplify the consideration we assume that the probability distribution p n does not depend on the distribution for 
the velocity. The simplest choice for the distribution of the number of particles in nearest surrounding is the Poisson 
distribution: 

Pn^^e- 1 , (22) 

TV. 

where X is mean number (density) of particles in the nearest surrounding. In order to get the analytic result in such a 
case it is expedient to use the non normalized order parameter instead of (19): 

u,= ^v 7 . (23) 

(hj) 

Then we get simple expression for the characteristic function Ga ■ 

G„(k) = e - MG ^. (24) 

The standard relation between the moments of the distribution and the derivatives of the characteristic function at k 
leads to: 

(u)=A(v)^(u) = (v). (25) 




FIGURE 3. Squares are experimental points, solid line is the solution Eq. (27). 

Note that such a simple relation between the order parameter and the particle velocity is due to specific form of the 
Poisson distribution. It allows to justify the expression for the relaxation rate Eq. (16) in low density limit. The velocity 
distribution function depends on the specific way of introducing the source of noise. Below we consider two types of 
noise which are widely used in simulations. 



3.1. Scalar noise 

Here we consider the case of so called scalar noise. This is the most obvious way to introduce the stochastic source 
into the dynamic of a particle. At every step the random increment of the angle is added. Using the results obtained 
above for the velocity distribution function, we choose it in accordance with Eq. (18) 

f y =Ce X % cosa , C = \— Y . (26) 

2kIq (A 5) 

u = )°{ (27) 

It's obvious, that it has trivial solution u = 0, which losses its stability depending on the average density X and the 
diffusion coefficient D (i.e. the intensity of scalar noise). Expanding Eq. (27) near the trivial solution we obtain : 



Q ={m- l ) u+ i^ u+o{u) - (28) 

From here the critical density value is as following: 

X c = 2D (29) 

The comparison of the solution of Eq. (27) and the results of numerical simulation obtained in [6, 17] is on Fig. 3. From 
Eq. (20) it clear that behavior of the order parameter near the critical value is determined by the analytical properties 
of the characteristic function G„. Near the critical threshold the dependence of the order parameter has typical for the 
mean-field approximation square root dependence: 

uoc (A t -A) 1/2 . (30) 

The differentiability properties of the characteristic function Eq. (20) and therefore the applicability of the expansion 
Eq. (28) depends on the character of the distribution of the neighbors. If the fluctuations of the density near the 
threshold value are big this can lead to slow convergence of p„. In such a case one can expect non smooth dependence 



of G u on the parameters of the distribution function / v , in particular on the average value of the local order parameter 
u. 

In general one need to construct the kinetic equation for the distribution function. Some attempts to derive such 
equation have been made in a way similar to classic Boltzmann approach [18] though only binary collisions were 
taken into account. This approximation is valid only for the system of low density. The applicability of these result to 
the systems of Vcsek's type is problematic because of the multiparticle character of the "collision" process. 



3.2. Vector noise 



The vectorial noise was introduced in [16] as another realistic model of the noisy environment. In such a case the 
random vector | is added, either to the local order parameter u [16] or u [19, 20]. Then the corresponding direction 
for the velocity of the particle is determined: 



0,(n + l)=Arg(u ( (n) + ^.; 



(31) 



In addition, the amplitude of the noise can be chosen so that <i; ; = £,Nj [16]. The results obtained in [16] revealed the 
difference between the VM with the scalar noise and raised the intensive discussion (see [11, 19, 20, 21]). 

Here we derive the one-particle velocity distribution function for the case of vectorial noise and show that it has 
essentially nonlinear character which leads to the apparent discontinuity in the dependence of the order parameter on 
the noise intensity. 

We assume that the distribution of the direction of the vector E, is uniform and independent on the the distribution 
of the number of neighbors, which is characterized by the corresponding probabilities p n . 
From simple geometrical consideration of vector noise algorithm it is easy to get the relation: 

1 (32) 
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If u < E, then the distribution function for the direction is as following: 
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(33) 



Thus the self-consistent equation for the order parameter is as following: 



u = { v)=F[ l 



(34) 



where 



= (cos a) = < 




if k/| < 1, 



if«/§ > 1. 



The solution of Eq. (34) is shown on Fig. 4. There is the interval of the noise intensity S,\ < % < §2, where two 
nontrivial solutions for the order parameter exist with the hysteretic (or subcritical) jump. It is clear that the branch 
where du/dE, > represents the unstable state in analogy with the situation typical for the first order phase transitions. 
For the model considered £i = 0.5 and I2 ~ 0.67. The situation here is analogous with that for the Kuramoto model 
[10]. In the latter case the type of bifurcation of the partially synchronized phase depends on the properties of the 
frequency distribution function g((0) of the oscillators, namely the sign of g"(0) (see [10]). Thus we state that the 
difference of the one-particle distribution function in case of scalar and vector noise is the source of the change the 
type of the bifurcation from supercritical for scalar noise to the subcritical for vector noise in the Vicsek model. This 
can explain the difference in the results of works [6] and [16]. 
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FIGURE 4. Solution of Eq. 34 



4. CONCLUSION 

In this paper we consider the continuum time limit for the Vicsek's algorithm [6] of the dynamics of the self propelling 
particles. It is shown that there is the time interval where the kinetic regime of the relaxation of the particle velocity to 
the local value of the average velocity of the neighbors takes place. The relaxation rate depends on the density linearly 
at least for not too big number of neighbors. The cases of vectorial and the scalar noises are considered. Within the 
proposed mean field model it is shown that the for the case of vectorial noise [16] the subcritical bifurcation of the 
stationary solution takes place. This is in contrast with the case of scalar noise originally considered in [6], where the 
supercritical transition occurs. 
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